---
title: "Moral suasion and the private provision of public goods: Evidence from the COVID-19 pandemic"
author: "Bos, B., Drupp, M.A., Meya, J.N., Quaas, M.F. (2020)"
output:
  html_document:
    df_print: paged
    code_folding: hide
    highlight: tango
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(fig.align = "center", warning = FALSE)

library(tidyverse)
library(magrittr)
library(here)
library(scales)
library(RColorBrewer)
library(ggpubr)
library(tikzDevice)

# Import data
survey_data <- read_csv(here::here("1_Data", "MoralSuasion_survey_data.csv"))

# Relabel and reorder the moral treatment variable
survey_data %<>% 
  mutate(moral_treatment_label = as.factor(case_when(
    moral_treatment == 1 ~ "Control",
    moral_treatment == 2 ~ "Deont.",
    moral_treatment == 3 ~ "Conseq.")))

survey_data$moral_treatment_label <- reorder(survey_data$moral_treatment_label, survey_data$moral_treatment)
levels(survey_data$moral_treatment_label)
table(survey_data$moral_treatment_label)

# Filter for the relevant respondents
survey_data %<>% 
  filter(time_required >= 3,
         time_required <= 60,
         valid == 1)

# Define the output path:
output_path = here::here("3_Graphs")
```


## Figure 2: Average planned private public good contributions by treatment group and treatment effects on the support for governmental regulation

```{r}
# Define a function for the plots
plot_means_by_treatment <- function(data, treatment_var, outcome, ymin, ymax, break_steps, ylab_label, label) {
  
  treatment_var <- enquo(treatment_var)
  outcome <- enquo(outcome)

  plot <- data %>% 
    group_by(!!treatment_var) %>% 
    summarise(mean = mean(!!outcome, na.rm = TRUE),
              sd = sd(!!outcome, na.rm = TRUE),
              se = sd / sqrt(n())) %>% 
    ggplot(aes(x = !!treatment_var, y = mean)) +
    geom_hline(yintercept = 0, color = "lightgray") +
    geom_bar(stat = "identity", width=.7, fill = "darkgrey") +
    geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width=.2)

  if (!missing(ymin) & !missing(ymax) & !missing(break_steps)) {
  plot <- plot +
    scale_y_continuous(limits = c(ymin, ymax),
                       breaks = seq(ymin, ymax, break_steps),
                       oob = rescale_none,
                       sec.axis = dup_axis(name = "", labels = NULL))
  }
  
  plot <- plot +
    ylab(ylab_label)+
    xlab("") + 
    theme_minimal() + 
    theme(
      # Adjust the grid lines
      panel.grid = element_blank(),
      # Increase the size of the axis and legend labels
      axis.text.x = element_text(size = 20, colour = "black", margin = unit(c(0.5,0.5,-0.5,0.5), "cm")),
      axis.text.y = element_text(size = 16, margin = unit(c(0.5,0.5,0.5,0.5), "cm")),
      axis.title.y = element_text(size = 20),
      # Add a border and axis ticks facing inwards
      panel.border = element_rect(colour = "black", fill=NA, size = 1.5),
      axis.ticks.y = element_line(),
      axis.ticks.length = unit(-0.1, "cm")
)
  
  return(plot)
}

# To test
# plot_means_by_treatment(data = survey_data,
#                         treatment_var = moral_treatment_label,
#                         outcome = change_contacts_next,
#                         ylab_label = "Relative to normal $(=1)$",
#                         ymin = 3, ymax = 4, break_steps = 0.5) +
#   stat_compare_means(data = survey_data, aes(x = moral_treatment_label, y = change_contacts_next, label = ..p.signif..),
#                      method = "t.test", ref.group = "Control",
#                      label.y = 4,
#                      symnum.args = list(cutpoints = c(0, 0.01, 0.05, 0.1, 1), symbols = c("***", "**", "*", "")))
```


```{r}
# a) Planned contacts
(plot1 <- plot_means_by_treatment(data = survey_data,
                        treatment_var = moral_treatment_label,
                        outcome = change_contacts_next,
                        ylab_label = "Relative to normal",
                        ymin = 3, ymax = 4, break_steps = 0.5) +
  stat_compare_means(data = survey_data, aes(x = moral_treatment_label, y = change_contacts_next, label = ..p.signif..),
                     method = "t.test", ref.group = "Control",
                     label.y = 4,
                     symnum.args = list(cutpoints = c(0, 0.01, 0.05, 0.1, 1), symbols = c("***", "**", "*", "")))
)


ggsave(filename = "fig_2_contacts.pdf", plot = plot1,
       path = output_path,
       width = 7, height = 5)

# Save for Latex
tikz(file = paste0(output_path, "/fig_2_contacts.tex"), width = 7, height = 5)
  plot1
dev.off()


# b) Planned hand cleaning effort
(plot2 <- plot_means_by_treatment(data = survey_data,
                        treatment_var = moral_treatment_label,
                        outcome = change_cleaning_hands_next,
                        ylab_label = "Relative to normal",
                        ymin = 10, ymax = 12, break_steps = 1) +
  stat_compare_means(data = survey_data, aes(x = moral_treatment_label, y = change_cleaning_hands_next, label = ..p.signif..),
                     method = "t.test", ref.group = "Control",
                     label.y = 12,
                     symnum.args = list(cutpoints = c(0, 0.01, 0.05, 0.1, 1), symbols = c("***", "**", "*", "")))
)

# Save as pdf
ggsave(filename = "fig_2_hc.pdf", plot = plot2,
       path = output_path,
       width = 7, height = 5)

# Save for Latex
tikz(file = paste0(output_path, "/fig_2_hc.tex"), width = 7, height = 5)
  plot2
dev.off()


# c) Support for gov. regulation
(plot3 <- survey_data %>% 
    group_by(moral_treatment_label) %>% 
    summarise(mean = mean(z_support, na.rm = TRUE),
              sd = sd(z_support, na.rm = TRUE),
              se = sd / sqrt(n())) %>% 
    ungroup() %>% 
    mutate(diff_to_mean = mean - first(mean)) %>% 
    ggplot(aes(x = moral_treatment_label, y = diff_to_mean)) +
    geom_hline(yintercept = 0, color = "lightgray") +
    geom_bar(stat = "identity", width=.7, fill = "darkgrey") +
    geom_errorbar(aes(ymin = diff_to_mean - se, ymax = diff_to_mean + se), width=.2) +
    ylab("Treatment effect") +
    xlab("") + 
    theme_minimal() + 
    theme(
      # Adjust the grid lines
      panel.grid = element_blank(),
      # Increase the size of the axis and legend labels
      axis.text.x = element_text(size = 20, colour = "black", margin = unit(c(0.5,0.5,-0.5,0.5), "cm")),
      axis.text.y = element_text(size = 16, margin = unit(c(0.5,0.5,0.5,0.5), "cm")),
      axis.title.y = element_text(size = 20),
      # Add a border and axis ticks facing inwards
      panel.border = element_rect(colour = "black", fill=NA, size = 1.5),
      axis.ticks.y = element_line(),
      axis.ticks.length = unit(-0.1, "cm")) +
  stat_compare_means(data = survey_data, aes(x = moral_treatment_label, y = z_support, label = ..p.signif..),
                     method = "t.test", ref.group = "Control",
                     label.y = 0.2,
                     symnum.args = list(cutpoints = c(0, 0.01, 0.05, 0.1, 1), symbols = c("***", "**", "*", "")))
)

ggsave(filename = "fig_2_support.pdf", plot = plot3,
       path = output_path,
       width = 7, height = 5)

# Save for Latex
tikz(file = paste0(output_path, "/fig_2_support.tex"), width = 7, height = 5)
  plot3
dev.off()


# d) Change in contacts wrt. gov. regulation
(plot4 <- survey_data %>% 
    group_by(moral_treatment_label) %>% 
    summarise(mean = mean(z_change_wrt_gov, na.rm = TRUE),
              sd = sd(z_change_wrt_gov, na.rm = TRUE),
              se = sd / sqrt(n())) %>% 
    ungroup() %>% 
    mutate(diff_to_mean = mean - first(mean)) %>% 
    ggplot(aes(x = moral_treatment_label, y = diff_to_mean)) +
    geom_hline(yintercept = 0, color = "lightgray") +
    geom_bar(stat = "identity", width=.7, fill = "darkgrey") +
    geom_errorbar(aes(ymin = diff_to_mean - se, ymax = diff_to_mean + se), width=.2) +
    ylab("Treatment effect") +
    xlab("") + 
    theme_minimal() + 
    theme(
      # Adjust the grid lines
      panel.grid = element_blank(),
      # Increase the size of the axis and legend labels
      axis.text.x = element_text(size = 20, colour = "black", margin = unit(c(0.5,0.5,-0.5,0.5), "cm")),
      axis.text.y = element_text(size = 16, margin = unit(c(0.5,0.5,0.5,0.5), "cm")),
      axis.title.y = element_text(size = 20),
      # Add a border and axis ticks facing inwards
      panel.border = element_rect(colour = "black", fill=NA, size = 1.5),
      axis.ticks.y = element_line(),
      axis.ticks.length = unit(-0.1, "cm")) +
  stat_compare_means(data = survey_data, aes(x = moral_treatment_label, y = z_support, label = ..p.signif..),
                     method = "t.test", ref.group = "Control",
                     label.y = 0.2,
                     symnum.args = list(cutpoints = c(0, 0.01, 0.05, 0.1, 1), symbols = c("***", "**", "*", "")))
)

# Save as pdf  
ggsave(filename = "fig_2_change_wrt_gov.pdf", plot = plot4,
       path = output_path,
       width = 7, height = 5)

# Save for Latex
tikz(file = paste0(output_path, "/fig_2_change_wrt_gov.tex"), width = 7, height = 5)
  plot4
dev.off()
```



## Figure A1: Distribution of outcome variables by treatment group
```{r}
# Define a function for these plots
plot_distribution <- function(data, treatment_var, outcome) {
  
  treatment_var <- enquo(treatment_var)
  outcome <- enquo(outcome)
  
  # Look up how many unique values there are
  num_bins <- data %>% 
    select(!!outcome) %>% 
    filter(!is.na(!!outcome)) %>%
    n_distinct()
  
  plot <- ggplot(data, aes(x = !!outcome,
                           y = after_stat(density),
                           fill = !!treatment_var)) +
  geom_histogram(bins = num_bins, position="dodge", center = 0) +
  scale_x_continuous(breaks = seq(1, num_bins, 1)) +
  scale_fill_manual(values = brewer.pal(3, "Paired")) +
  scale_y_continuous(expand = c(0,0),
                     sec.axis = dup_axis(name = "", labels = NULL)) +  
  ylab("Density") +
  xlab("") +
  theme_minimal() +
  theme(
    # Adjust the grid lines
    panel.grid = element_blank(),
    # Adjust the legend
    legend.position = "bottom",
    legend.title = element_blank(),
    legend.spacing.x = unit(0.5, 'cm'),
    # Increase the size of the axis and legend labels
    axis.text.x = element_text(size = 18, margin=unit(c(0.5,0.5,-0.75,0.5), "cm")),
    axis.text.y = element_text(size = 18, margin=unit(c(0.5,0.5,0.5,0.5), "cm")),
    axis.title.y = element_text(size = 20),
    legend.text = element_text(size = 20),
    # Add a border and axis ticks facing inwards
    panel.border = element_rect(colour = "black", fill=NA),
    axis.ticks.y = element_line(),
    axis.ticks.length = unit(-0.1, "cm")
    )

  return(plot)
}

# To test:
# plot_distribution(data = survey_data,
#                   treatment_var = moral_treatment_label,
#                   outcome = change_contacts_next)
```


```{r warning=FALSE}
# Note, graphs do not include rows with missing values in those variables.
# Hence, we turn off those warnings.

# a) Planned contacts
fig_distribution_contacts <- plot_distribution(data = survey_data,
                                               treatment_var = moral_treatment_label,
                                               outcome = change_contacts_next)

fig_distribution_contacts +
  labs(title = "a) Planned contacts")
  

ggsave(fig_distribution_contacts,
       filename = "fig_A1_contacts.pdf",
       path = output_path,
       width = 7, height = 5)

# Save for Latex
tikz(file = paste0(output_path, "/fig_A1_contacts.tex"), width = 7, height = 5)
  fig_distribution_contacts
dev.off()


# b) Planned hand cleaning effort
fig_distribution_hc <- plot_distribution(data = survey_data,
                                         treatment_var = moral_treatment_label,
                                         outcome = change_cleaning_hands_next)

fig_distribution_hc +
  labs(title = "b) Planned hand cleaning effort")


ggsave(fig_distribution_hc,
       filename = "fig_A1_hc.pdf",
       path = output_path,
       width = 7, height = 5)

# Save for Latex
tikz(file = paste0(output_path, "/fig_A1_hc.tex"), width = 7, height = 5)
  fig_distribution_hc
dev.off()


# c) Support for gov. regulation
fig_distribution_support <- plot_distribution(data = survey_data,
                                              treatment_var = moral_treatment_label,
                                              outcome = perception_gov_reg)

fig_distribution_support +
  labs(title = "c) Support for gov. regulation")


ggsave(fig_distribution_support,
       filename = "fig_A1_support.pdf",
       path = output_path,
       width = 7, height = 5)

# Save for Latex
tikz(file = paste0(output_path, "/fig_A1_support.tex"), width = 7, height = 5)
  fig_distribution_support
dev.off()



# d) Change in contacts wrt. gov. regulation
fig_distribution_change_wrt_gov <- plot_distribution(data = survey_data,
                                                     treatment_var = moral_treatment_label,
                                                     outcome = change_con_wrt_gov_reg)

fig_distribution_change_wrt_gov +
  labs(title = "d) Change in contacts wrt. gov. regulation")


ggsave(fig_distribution_change_wrt_gov,
       filename = "fig_A1_change_wrt_gov.pdf",
       path = output_path,
       width = 7, height = 5)

# Save for Latex
tikz(file = paste0(output_path, "/fig_A1_change_wrt_gov.tex"), width = 7, height = 5)
  fig_distribution_change_wrt_gov
dev.off()
```

